C $TITLE: 'FACTRS'
C $NOFLOATCALLS
      SUBROUTINE FACTRS(A,SCRATC,NP,NROW,IP,IX,IU1,IU2,IU3,IU4,
     1 LD2,IRESRV)
C
C     FACTRS, FOR SYMMETRIC STRUCTURE, TRANSFORMS SUBMATRICIES TO FORM
C     MATRICIES OF THE SYMMETRIC MODES AND CALLS ROUTINE TO FACTOR
C     MATRICIES.  IF NO SYMMETRY, THE ROUTINE IS CALLED TO FACTOR THE
C     COMPLETE MATRIX.
C
      INTEGER*4 I2,NP,NROW,IDM1,KA,KK,NOP
      INTEGER*4 IMAT,NPBLK,NLAST,NLSYM,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
CLARGE: A
      COMPLEX A
      COMPLEX*16 SCRATC
      COMMON/MATPAR/ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,IMAT,
     1 ICASX,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
      DIMENSION A(1),SCRATC(LD2),IP(LD2),IX(LD2)
C**
C     D     WRITE(*,*) '  FACTRS: ICASE=',ICASE,' ICASX=',ICASX
      IRESRV=IRESRV
C**
      IDM1=1
      NOP=NROW/NP
      IF (ICASE.GT.2) GO TO 2
C**
      DO 1 KK=1,NOP
      KA=(KK-1)*NP+1
C**
C     D     WRITE(*,*) '  FACTRS: CALL FACTR, KA=',KA,' NP=',NP,' NROW=',NROW
C**
      CALL FACTR(A(KA),SCRATC,NP,NROW,IP(KA),LD2)
C**
C     D     WRITE(*,*) '  FACTRS: RTRN FACTR'
C**
1     CONTINUE
C**
C     D     WRITE(*,*) '  FACTRS: RETURN AFTER 1'
C**
      RETURN
2     IF (ICASE.GT.3) GO TO 3
C
C     FACTOR SUBMATRICIES, OR FACTOR COMPLETE MATRIX IF NO SYMMETRY
C     EXISTS.
C**
C     D     WRITE(*,*) '  FACTRS: CALL FACIO'
C**
      CALL FACIO (A,SCRATC,NROW,NOP,IX,IU1,IU2,IU3,IU4,LD2)
C**
C     D     WRITE(*,*) '  FACTRS: RTRN FACIO'
C     D     WRITE(*,*) '  FACTRS: CALL LUNSCR'
C**
      CALL LUNSCR (A,NROW,NOP,IP,IX,IU2,IU3,IU4)
C**
C     D     WRITE(*,*) '  FACTRS: RTRN LUNSCR'
C     D     WRITE(*,*) '  FACTRS: RETURN BEFORE 3'
C**
      RETURN
C
C     REWRITE THE MATRICES BY COLUMNS ON TAPE 13
C
3     I2=2*NPBLK*NROW
      REWIND IU2
      DO 5 K=1,NOP
      REWIND IU1
      ICOLS=NPBLK
      IR2=K*NP
      IR1=IR2-NP+1
      DO 5 L=1,NBLOKS
      IF (NBLOKS.EQ.1.AND.K.GT.1) GO TO 4
C**
C     D     WRITE(*,*) '  FACTRS: CALL BLCKIN'
C**
      CALL BLCKIN (A,IDM1,I2,1,602,IU1)
C**
C     D     WRITE(*,*) '  FACTRS: RTRN BLCKIN'
C**
      IF (L.EQ.NBLOKS) ICOLS=NLAST
4     IRR1=IR1
      IRR2=IR2
      DO 5 ICOLDX=1,ICOLS
      WRITE (IU2) (A(I),I=IRR1,IRR2)
      IRR1=IRR1+NROW
      IRR2=IRR2+NROW
5     CONTINUE
      REWIND IU1
      REWIND IU2
      IF (ICASE.EQ.5) GO TO 8
      REWIND IU3
      IRR1=NP*NP
      DO 7 KK=1,NOP
      IR1=1-NP
      IR2=0
      DO 6 I=1,NP
      IR1=IR1+NP
      IR2=IR2+NP
6     READ (IU2) (A(J),J=IR1,IR2)
      KA=(KK-1)*NP+1
C**
C     D     WRITE(*,*) '  FACTRS: CALL FACTR AFTER 6, KA=',KA,' NP=',NP
C**
      CALL FACTR(A,SCRATC,NP,NP,IP(KA),LD2)
C**
C     D     WRITE(*,*) '  FACTRS: RTRN FACTR AFTER 6'
C**
      WRITE (IU3) (A(I),I=1,IRR1)
7     CONTINUE
      REWIND IU2
      REWIND IU3
      RETURN
8     I2=2*NPSYM*NP
      DO 10 KK=1,NOP
      J2=NPSYM
      DO 10 L=1,NBLSYM
      IF (L.EQ.NBLSYM) J2=NLSYM
      IR1=1-NP
      IR2=0
      DO 9 J=1,J2
      IR1=IR1+NP
      IR2=IR2+NP
9     READ (IU2) (A(I),I=IR1,IR2)
10     CONTINUE
C**
C     D     WRITE(*,*) '  FACTRS: CALL BLCKOT'
C**
      CALL BLCKOT (A,IDM1,I2,1,193,IU1)
C     D     WRITE(*,*) '  FACTRS: RTRN BLCKOT'
C**
      REWIND IU1
C**
C     D     WRITE(*,*) '  FACTRS: CALL FACIO'
C**
      CALL FACIO (A,SCRATC,NP,NOP,IX,IU1,IU2,IU3,IU4,LD2)
C**
C     D     WRITE(*,*) '  FACTRS: RTRN FACIO'
C     D     WRITE(*,*) '  FACTRS: CALL LUNSCR'
C**
      CALL LUNSCR (A,NP,NOP,IP,IX,IU2,IU3,IU4)
C**
C     D     WRITE(*,*) '  FACTRS: RTRN LUNSCR'
C     D     WRITE(*,*) '  FACTRS: RETURN AT END'
C**
      RETURN
      END
C
C
C
      SUBROUTINE REBLK (B,BX,NB,NBX,N2C)
C     REBLOCK ARRAY B IN N.G.F. SOLUTION FROM BLOCKS OF ROWS ON TAPE14
C     TO BLOCKS OF COLUMNS ON TAPE16
      INTEGER*4 NB,NBX,N2C
      INTEGER*4 IMAT,NPBLK,NLAST,NLSYM,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
CLARGE: B,BX
      COMPLEX B,BX
      DIMENSION B(NB,1),BX(NBX,1)
      COMMON/MATPAR/ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,IMAT,
     1 ICASX,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
C**
C     E     WRITE(*,*) '   REBLK: START NBBL=',NBBL
C**
      REWIND 16
      NIB=0
      NPB=NPBL
      DO 3 IB=1,NBBL
      IF (IB.EQ.NBBL) NPB=NLBL
      REWIND 14
      NIX=0
      NPX=NPBX
      DO 2 IBX=1,NBBX
      IF (IBX.EQ.NBBX) NPX=NLBX
      READ (14) ((BX(I,J),I=1,NPX),J=1,N2C)
      DO 1 I=1,NPX
      IX=I+NIX
      DO 1 J=1,NPB
1     B(IX,J)=BX(I,J+NIB)
2     NIX=NIX+NPBX
      WRITE (16) ((B(I,J),I=1,NB),J=1,NPB)
3     NIB=NIB+NPBL
      REWIND 14
      REWIND 16
C**
C     E     WRITE(*,*) '   REBLK: RETURN'
C**
      RETURN
      END
C
C
C
      SUBROUTINE FACIO (A,D,NROW,NOP,IP,IU1,IU2,IU3,IU4,LD2)
C
C     FACIO CONTROLS I/O FOR OUT-OF-CORE FACTORIZATION
C
      INTEGER*4 I1,I2,I3,I4,IT,NROW
CLARGE A
      COMPLEX A
      COMPLEX*16 D
      DIMENSION A(NROW,1),IP(NROW),D(LD2)
      INTEGER*4 NOP,KK
      INTEGER*4 IMAT,NPBLK,NLAST,NLSYM,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
      COMMON/MATPAR/ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,IMAT,
     1 ICASX,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
C**
C     D     WRITE(*,*) '   FACIO: START'
C**
      IT=2*NPSYM*NROW
      NBM=NBLSYM-1
      I1=1
      I2=IT
      I3=I2+1
      I4=2*IT
      TIME=0.
      REWIND IU1
      REWIND IU2
      DO 3 KK=1,NOP
      KA=(KK-1)*NROW+1
      IFILE3=IU1
      IFILE4=IU3
      DO 2 IXBLK1=1,NBM
      REWIND IU3
      REWIND IU4
      CALL BLCKIN (A,I1,I2,1,17,IFILE3)
      IXBP=IXBLK1+1
      DO 1 IXBLK2=IXBP,NBLSYM
      CALL BLCKIN (A,I3,I4,1,18,IFILE3)
      CALL SECOND (T1)
      CALL LFACTR (A,D,NROW,IXBLK1,IXBLK2,IP(KA),LD2)
      CALL SECOND (T2)
      TIME=TIME+T2-T1
      IF (IXBLK2.EQ.IXBP) CALL BLCKOT (A,I1,I2,1,19,IU2)
      IF (IXBLK1.EQ.NBM.AND.IXBLK2.EQ.NBLSYM) IFILE4=IU2
      CALL BLCKOT (A,I3,I4,1,20,IFILE4)
1     CONTINUE
      IFILE3=IU3
      IFILE4=IU4
      IF ((IXBLK1/2)*2.NE.IXBLK1) GO TO 2
      IFILE3=IU4
      IFILE4=IU3
2     CONTINUE
3     CONTINUE
      REWIND IU1
      REWIND IU2
      REWIND IU3
      REWIND IU4
      WRITE(*,4)  TIME
      RETURN
C
4     FORMAT (' CP TIME TAKEN FOR FACTORIZATION = ',1P,E12.5)
      END
C
C
C
      SUBROUTINE LUNSCR (A,NROW,NOP,IX,IP,IU2,IU3,IU4)
C
C     S/R WHICH UNSCRAMBLES, SCRAMBLED FACTORED MATRIX
C
      INTEGER*4 I1,I2,NROW
CLARGE A
      COMPLEX A
      COMPLEX*16 TEMP
      DIMENSION A(NROW,1),IP(NROW),IX(NROW)
      INTEGER*4 NOP,KA,KK
      INTEGER*4 IMAT,NPBLK,NLAST,NLSYM,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
      COMMON/MATPAR/ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,IMAT,
     1 ICASX,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
C**
C     D     WRITE(*,*) '   LUNSCR: START'
C**
      I1=1
      I2=2*NPSYM*NROW
      NM1=NROW-1
      REWIND IU2
      REWIND IU3
      REWIND IU4
      DO 9 KK=1,NOP
      KA=(KK-1)*NROW
      DO 4 IXBLK1=1,NBLSYM
      CALL BLCKIN (A,I1,I2,1,121,IU2)
      K1=(IXBLK1-1)*NPSYM+2
      IF (NM1.LT.K1) GO TO 3
      J2=0
      DO 2 K=K1,NM1
      IF (J2.LT.NPSYM) J2=J2+1
      IPK=IP(K+KA)
      DO 1 J=1,J2
      TEMP=A(K,J)
      A(K,J)=A(IPK,J)
      A(IPK,J)=TEMP
1     CONTINUE
2     CONTINUE
3     CONTINUE
      CALL BLCKOT (A,I1,I2,1,122,IU3)
4     CONTINUE
      DO 5 IXBLK1=1,NBLSYM
      BACKSPACE IU3
      IF (IXBLK1.NE.1) BACKSPACE IU3
      CALL BLCKIN (A,I1,I2,1,123,IU3)
      CALL BLCKOT (A,I1,I2,1,124,IU4)
5     CONTINUE
      DO 6 I=1,NROW
      IX(I+KA)=I
6     CONTINUE
      DO 7 I=1,NROW
      IPI=IP(I+KA)
      IXT=IX(I+KA)
      IX(I+KA)=IX(IPI+KA)
      IX(IPI+KA)=IXT
7     CONTINUE
      IF (NOP.EQ.1) GO TO 9
      NB1=NBLSYM-1
C     SKIP NB1 LOGICAL RECORDS FORWARD
      DO 8 IXBLK1=1,NB1
      CALL BLCKIN (A,I1,I2,1,125,IU3)
8     CONTINUE
9     CONTINUE
      REWIND IU2
      REWIND IU3
      REWIND IU4
      RETURN
      END
C
C
C
      SUBROUTINE LFACTR (A,D,NROW,IX1,IX2,IP,LD2)
C
C     LFACTR PERFORMS GAUSS-DOOLITTLE MANIPULATIONS ON THE TWO BLOCKS OF
C     THE TRANSPOSED MATRIX IN CORE STORAGE.  THE GAUSS-DOOLITTLE
C     ALGORITHM IS PRESENTED ON PAGES 411-416 OF A. RALSTON -- A FIRST
C     COURSE IN NUMERICAL ANALYSIS.  COMMENTS BELOW REFER TO COMMENTS IN
C     RALSTONS TEXT.
C
CLARGE A
      COMPLEX A
      COMPLEX*16 D,AJR
      REAL*8 DMAX,ELMAG
      INTEGER*4 NROW
      INTEGER R,R1,R2,PJ,PR
      LOGICAL L1,L2,L3
      INTEGER*4 IMAT,NPBLK,NLAST,NLSYM,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
      COMMON/MATPAR/ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,IMAT,
     1 ICASX,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
      DIMENSION A(NROW,1),IP(NROW),D(LD2)
C**
C     D     WRITE(*,*) '   LFACTR: START'
C**
      IFLG=0
C
C     INITIALIZE R1,R2,J1,J2
C
      L1=IX1.EQ.1.AND.IX2.EQ.2
      L2=(IX2-1).EQ.IX1
      L3=IX2.EQ.NBLSYM
      IF (L1) GO TO 1
      GO TO 2
1     R1=1
      R2=2*NPSYM
      J1=1
      J2=-1
      GO TO 5
2     R1=NPSYM+1
      R2=2*NPSYM
      J1=(IX1-1)*NPSYM+1
      IF (L2) GO TO 3
      GO TO 4
3     J2=J1+NPSYM-2
      GO TO 5
4     J2=J1+NPSYM-1
5     IF (L3) R2=NPSYM+NLSYM
      DO 16 R=R1,R2
C
C     STEP 1
C
      DO 6 K=J1,NROW
      D(K)=A(K,R)
6     CONTINUE
C
C     STEPS 2 AND 3
C
      IF (L1.OR.L2) J2=J2+1
      IF (J1.GT.J2) GO TO 9
      IXJ=0
      DO 8 J=J1,J2
      IXJ=IXJ+1
      PJ=IP(J)
      AJR=D(PJ)
      A(J,R)=AJR
      D(PJ)=D(J)
      JP1=J+1
      DO 7 I=JP1,NROW
      D(I)=D(I)-A(I,IXJ)*AJR
7     CONTINUE
8     CONTINUE
9     CONTINUE
C
C     STEP 4
C
      J2P1=J2+1
      IF (L1.OR.L2) GO TO 11
      IF (NROW.LT.J2P1) GO TO 16
      DO 10 I=J2P1,NROW
      A(I,R)=D(I)
10    CONTINUE
      GO TO 16
11    DMAX=DREAL(D(J2P1)*DCONJG(D(J2P1)))
      IP(J2P1)=J2P1
      J2P2=J2+2
      IF (J2P2.GT.NROW) GO TO 13
      DO 12 I=J2P2,NROW
      ELMAG=DREAL(D(I)*DCONJG(D(I)))
      IF (ELMAG.LT.DMAX) GO TO 12
      DMAX=ELMAG
      IP(J2P1)=I
12    CONTINUE
13    CONTINUE
      IF (DMAX.LT.1.E-10) IFLG=1
      PR=IP(J2P1)
      A(J2P1,R)=D(PR)
      D(PR)=D(J2P1)
C
C     STEP 5
C
      IF (J2P2.GT.NROW) GO TO 15
      AJR=1.D0/A(J2P1,R)
      DO 14 I=J2P2,NROW
      A(I,R)=D(I)*AJR
14    CONTINUE
15    CONTINUE
      IF (IFLG.EQ.0) GO TO 16
      WRITE(*,17)  J2,DMAX
      IFLG=0
16    CONTINUE
C**
C     D     WRITE(*,*) '   LFACTR: RETURN'
C**
      RETURN
C
17    FORMAT (1H ,6HPIVOT(,I3,2H)=,1P,D16.8)
      END
C
C
C
      SUBROUTINE FACGF(A,B,C,D,BX,SCRATC,IP,IX,NP,N1,
     1 MP,M1,N1C,N2C,LD2,IRESRV)
C     FACGF COMPUTES AND FACTORS D-C(INV(A)B).
      INTEGER*4 NP,N1,MP,M1,N1C,N2C,IDM2,N1CP
      INTEGER*4 IMAT,NPBLK,NLAST,NLSYM,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
CLARGE: A,B,C,D,BX
      COMPLEX A,B,C,D,BX
      COMPLEX*16 SUM,SCRATC
      COMMON/MATPAR/ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,IMAT,
     1 ICASX,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
      DIMENSION A(1),B(N1C,1),C(N1C,1),D(N2C,1),BX(N1C,1),
     1 SCRATC(LD2),IP(LD2),IX(LD2)
      IF (N2C.EQ.0) RETURN
C**
C     D     WRITE(*,*) '  FACGF: START N1C=',N1C,' N2C=',N2C,' NBBL=',NBBL,
C     D    1 ' ICASE=',ICASE,' ICASX=',ICASX
C**
      IBFL=14
      IF (ICASX.LT.3) GO TO 1
C     CONVERT B FROM BLOCKS OF ROWS ON T14 TO BLOCKS OF COL. ON T16
C**
C     D     WRITE(*,*) '  FACGF: CALL REBLK'
C**
      CALL REBLK (B,C,N1C,NPBX,N2C)
C**
C     D     WRITE(*,*) '  FACGF: RTRN REBLK'
C**
      IBFL=16
1     NPB=NPBL
      IF (ICASX.EQ.2) REWIND 14
C     COMPUTE INV(A)B AND WRITE ON TAPE14
      DO 2 IB=1,NBBL
      IF (IB.EQ.NBBL) NPB=NLBL
      IF (ICASX.GT.1) READ (IBFL) ((BX(I,J),I=1,N1C),J=1,NPB)
C**
C     D     WRITE(*,*) '  FACGF: CALL SOLVES'
C**
      CALL SOLVES(A,BX,SCRATC,N1C,NP,N1,MP,M1,IP,NPB,13,13,LD2,
     1 IRESRV)
C**
C     D     WRITE(*,*) '  FACGF: RTRN SOLVES N1C=',N1C,' N2C=',N2C,
C     D    1' NBBL=',NBBL
C**
      IF (ICASX.EQ.2) REWIND 14
      IF (ICASX.GT.1) WRITE (14) ((BX(I,J),I=1,N1C),J=1,NPB)
2     CONTINUE
      IF (ICASX.EQ.1) GO TO 3
      REWIND 11
      REWIND 12
      REWIND 15
      REWIND IBFL
3     NPC=NPBL
C     COMPUTE D-C(INV(A)B) AND WRITE ON TAPE11
      DO 8 IC=1,NBBL
      IF (IC.EQ.NBBL) NPC=NLBL
      IF (ICASX.EQ.1) GO TO 4
      READ (15) ((C(I,J),I=1,N1C),J=1,NPC)
      READ (12) ((D(I,J),I=1,N2C),J=1,NPC)
      REWIND 14
4     NPB=NPBL
      NIC=0
      DO 7 IB=1,NBBL
      IF (IB.EQ.NBBL) NPB=NLBL
      IF (ICASX.GT.1) READ (14) ((B(I,J),I=1,N1C),J=1,NPB)
      DO 6 I=1,NPB
      II=I+NIC
      DO 6 J=1,NPC
      SUM= CMPLX(0.,0.)
      DO 5 K=1,N1C
5     SUM=SUM+B(K,I)*C(K,J)
6     D(II,J)=D(II,J)-SUM
7     NIC=NIC+NPBL
      IF (ICASX.GT.1) WRITE (11) ((D(I,J),I=1,N2C),J=1,NPBL)
8     CONTINUE
      IF (ICASX.EQ.1) GO TO 9
      REWIND 11
      REWIND 12
      REWIND 14
      REWIND 15
9     N1CP=N1C+1
C     FACTOR D-C(INV(A)B)
      IF (ICASX.GT.1) GO TO 10
C**
C     D     WRITE(*,*) '  FACGF: CALL FACTR'
C**
      CALL FACTR(D,SCRATC,N2C,N2C,IP(N1CP),LD2)
C**
C     D     WRITE(*,*) '  FACGF: RTRN FACTR N1C=',N1C,' N2C=',N2C,
C     D    1' NBBL=',NBBL
C**
      GO TO 13
10    IF (ICASX.EQ.4) GO TO 12
      NPB=NPBL
      IC=0
C**
      DO 11 IB=1,NBBL
      IF (IB.EQ.NBBL) NPB=NLBL
      II=IC+1
      IC=IC+N2C*NPB
C**
C     D       WRITE(*,*) ' II=',II,' IB=',IB,' IC=',IC,' NPB=',NPB
C     D       WRITE(*,*) ' (B(I,1),I=II,IC)=',(B(I,1),I=II,IC)
C**
11    READ (11) (B(I,1),I=II,IC)
      REWIND 11
C**
C     D     WRITE(*,*) '  FACGF: CALL FACTR AFTER 11'
C**
      CALL FACTR(B,SCRATC,N2C,N2C,IP(N1CP),LD2)
C**
C     D     WRITE(*,*) '  FACGF: RTRN FACTR AFTER 11'
C**
      NIC=N2C*N2C
      WRITE (11) (B(I,1),I=1,NIC)
      REWIND 11
      GO TO 13
C**
12    NBLSYS=NBLSYM
      NPSYS=NPSYM
      NLSYS=NLSYM
      ICASS=ICASE
      NBLSYM=NBBL
      NPSYM=NPBL
      NLSYM=NLBL
      ICASE=3
      IDM2=1
C**
C     D     WRITE(*,*) '  FACGF: CALL FACIO'
C**
      CALL FACIO (B,SCRATC,N2C,IDM2,IX(N1CP),11,12,16,11,LD2)
C**
C     D     WRITE(*,*) '  FACGF: RTRN FACIO'
C     D     WRITE(*,*) '  FACGF: CALL LUNSCR'
C**
      CALL LUNSCR (B,N2C,IDM2,IP(N1CP),IX(N1CP),12,11,16)
C**
C     D     WRITE(*,*) '  FACGF: RTRN LUNSCR'
C**
      NBLSYM=NBLSYS
      NPSYM=NPSYS
      NLSYM=NLSYS
      ICASE=ICASS
13     CONTINUE
C**
C     D     WRITE(*,*) '  FACGF: RETURN AT END ICASE=',ICASE
C**
      RETURN
      END
